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We analyse the angular dynamics of a neutrally buoyant nearly spherical particle immersed in a 
steady general linear flow. The hydrodynamic torque acting on the particle is obtained by means 
of a reciprocal theorem, regular perturbation theory exploiting the small eccentricity of the nearly 
spherical particle, and assuming that inertial effects are small, but finite. 

PACS numbers: 83.10.Pp,47.15.G-,47.55.Kf,47.10.-g 


I. INTRODUCTION 

In this article we derive an effective equation of mo¬ 
tion for the orientational dynamics of a neutrally buoy¬ 
ant, nearly spherical axisymmetric particle suspended in 
a time-independent linear flow. Our result is valid to 
leading order in the shear Reynolds number Re^ and the 
particle eccentricity e. Terms of order Re^, e^, and e^Re^ 
are neglected. 

Our motivation was two-fold. Firstly, we have recently 
computed the stability of log-rolling and tumbling orbits 
of a neutrally buoyant spheroid in a simple shear at weak 
fluid and particle inertia mm- The calculations leading 
to these results are quite involved. We therefore decided 
to check our calculations by an alternative method, sum¬ 
marised below. We refer to Ref. [5] for a summary of 
the background of the problem, and for a discussion of 
the implications of the results. Secondly, the results de¬ 
scribed in the present article are valid for general linear 
flows while those given in Refs. mm pertain to the par¬ 
ticular (and important) case of a simple shear flow. We 
believe that it is of interest to obtain results for general 
linear flows because this can be a first step towards de¬ 
scribing the effect of a time-dependent but slowly varying 
perturbation of the flow. 


II. FORMULATION OF THE PROBLEM 



FIG. 1: Color online. Left: the surface of a nearly spherical 
particle is constructed as the surface of revolution of the curve 
Eq. Q. Right: simple shear. The flow-shear plane is spanned 
by ei and e 2 , and the vorticity points in the negative 63 - 
direction. The orientation of the vector ti pointing along 
the symmetry axis of the particle is expressed in spherical 
coordinates (0,4)). 


measures the effect of fluid inertia, the Stokes number 
St = (/9p/pf)Res measuring the importance of particle 
inertia, and the particle eccentricity e. Here pf and Pp 
denote, respectively, the densities of the fluid and of the 
particle, s denotes the shear rate of the linear flow, and 
p the dynamic viscosity of the fluid. 

Using the inverse shear rate as the time scale, the semi¬ 
axis length a as the length scale and ps as the pressure 
scale, the equations governing the angular dynamics of a 
neutrally buoyant particle read in the laboratory frame 
of reference 


We consider a nearly spherical particle corresponding 
to an ellipsoid of revolution (around the axis ti) of low 
eccentricity as depicted in Fig. The surface of the 
particle is parametrised as 

r{f)) = 1 - esin2(i?) - ^ sin2(^) cos'^[^) + 0{e^) . (1) 

Here f) denotes the polar angle made by any vector r with 
the orientation of the axis of revolution ti. The eccen¬ 
tricity e is a small parameter. Lengths are normalized by 
the semi-axis length a along the direction ti (Fig. [^. 

Suppose that the particle is immersed in a steady gen¬ 
eral linear flow. In this case the angular dynamics of 
the particle is governed by three dimensionless parame¬ 
ters: the shear Reynolds number Reg = sa^pf/p which 


dti , , 

( 2 ) 

and 

St^l^J • =Th. (3) 

In Eq. ([^, J denotes the moment-of-inertia tensor of 
the particle, lj its angular velocity, and Th is the hydro- 
dynamic torque acting on the particle. Throughout this 
paper the centred dot (•) defines a simply contracted ten¬ 
sor product, uppercase letters are used to denote matrix- 
tensors and lowercase letters are used for simple vectors. 

For a neutrally buoyant particle we have St = Reg. 
Making this substitution only at the end of the calcula- 
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tion allows us to separate the effects of particle and fluid 
inertia. 

The hydrodynamic torque is given by 



with r = X — Xp. 

We use the reciprocal theorem to determine the torque 
acting on the spheroid. The reciprocal theorem relates 
integrals of the velocity and stress fields of the two in¬ 
compressible Newtonian fluids: 


The integral is over the particle surface =5^, and n is the 
outward surface normal. For an incompressible Newto¬ 
nian fluid the stress tensor reads S = —p I -I- 2 S , where p 
is the pressure, I the identity tensor and S the symmetric 
part of the fluid-velocity gradient tensor. To determine 
the hydrodynamic torque acting on the spheroid, we use 
a reciprocal theorem mz]. 


III. METHOD 

We consider a general steady linear ambient flow. In 
the laboratory frame, it takes the form 


I (^(S-n)-u-(S-n)-SjdS' (5) 

5-(V-S)-u-(V-S))dF, 

where the integral on the r.h.s is over the entire volume 'f' 
outside the particle. In Eq. © the set (u, S) represents 
the solutions of the problem of interest and (u, S) are the 
solutions of a suitable auxiliary problem describing the 
creeping flow produced by a particle moving with angular 
velocity a; in an otherwise quiescent fluid. 

We also introduce the decompositions 

u = A°° • r + u' and S = -f S' , 



= A°° • X , 

where A°° a constant tensor. We decompose the gradient 
matrix A°° into its symmetric (S°°) and antisymmetric 
(0~) parts: 

A°° = -k . 

The anti-symmetric tensor is linked to the vorticity 
of the unperturbed flow through 

0°° • X = X X where i V x . 
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The vectors Xp and Up denote the the position and the 
velocity of the centre-of-mass of the particle. As the 
spheroid is assumed to be neutrally buoyant, we consider 
it to be advected along streamlines with a vanishing slip 
velocity 

Up = A°° • Xp . 


where u' and S' correspond to the perturbation flow and 
the perturbation stress tensor induced by the spheroid. 
Using these decompositions it follows from the bound¬ 
ary conditions for u' both on the particle surface and at 
infinity that: 

i^s-Th- (S°° • r) d^ (6) 

= a; • Th -I- Rcs / f(u')-udU. 

Jv 

In this equation we have introduced the slip angular ve¬ 
locity 

LJs = ^ — 1 

and the function f(u') in Eq. ^ is given by 
f(u') = ^ -p A~ • u' -P ((A'^ • r) • v)u' + u' • Vu'. 


We consider the equations governing the fluid motion 
in a translating frame of reference centred on the par¬ 
ticle. This frame being non-inertial, the pseudo-force 
—Re^dUp/dt appears in the force balance acting on the 
particle. However, it does not appear in the equations of 
the perturbation flow below. In the moving frame of ref¬ 
erence the Navier-Stokes equations for the fluid velocity 
u read: 

V-u = 0, 

Res (^ + (u-V)u) = V-S-Re,^, (4) 

with boundary conditions 


To derive Eq. § we have also used the fact that the 
torque due to the unperturbed stress vanishes. To see 
this recall that in the laboratory frame of reference, the 
unperturbed flow field satisfies the Navier-Stokes equa¬ 
tions A^ • X = V • Using Stokes integration theorem 
it follows that the torque due to the unperturbed stress 
tensor is given by 



( 7 ) 


where the integral is over the volume 'fp inside the par¬ 
ticle. The symmetrical shape of the particle implies that 
the first integral on the r.h.s of Eq. 0 is zero. Incom¬ 
pressibility and the vorticity equation for a steady flow 


u = ci;xr, r G ^ and u — )■ A°° • r , r —)■ oo , 


Tr(S~) =0 and = 0 


( 8 ) 
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imply that the second integral also vanishes. 

Eq. (§ allows to compute the hydrodynamic torque 
Tfi provided that we can determine the contribution of 
the integral over the entire volume f' outside the parti¬ 
cle. To leading order HE] one can replace the actual 
perturbation flow u' in the integral by its limit of vanish¬ 
ingly small ROg. In other words, u' can be replaced by 
the solution of 


In order to determine the second term on the r.h.s of 
Eq. Q, which is, as already said, the most difficult part 
of this investigation, we first make use of the fact that 
the slip angular velocity is expected to be of order 0(e). 
Indeed, Jeffery’s theory El predicts 

lim = Ati X (S°° • ti) . 

Re,^o 


0 = V • u', 0 = -Vp' + Au' (9) 


with boundary conditions 


u' = (o) — a;°°) X r — • r for r S , 

u'—> 0 as r —>• oo . (10) 


To determine u' we follow Refs. |S1[5] and use a perturba¬ 
tion method that assumes that the particle eccentricity e 
is small. Our approach is based on the general Lamb so¬ 
lution [To] of the Stokes problem ( 9|10 | derived in terms 
of spherical harmonics. Details are given in the appendix. 


IV. RESULTS FOR GENERAL LINEAR FLOWS 


Eq. @ makes it possible to derive the torque acting 
on a neutrally buoyant spheroid immersed in a general 
linear flow. For the auxiliary torque - the first term on 
the l.h.s. of Eq. (§ - we find: 

„ 9 459 n\.^ 

The second term on the l.h.s. of Eq. (§ evaluates to: 


S • n) • (S°° • r) d5 

= 8^(e-^e2)((S-.ti)xti) 


•w + 0(e3). 


The factor A is a shape factor determined by the particle 
aspect ratio A (the ratio between the length of the particle 
along its symmetry axis and its length transverse the 
symmetry axis) through 


_ 1 


Using the fact that the slip angular velocity is of order 
0{e) and making use of the identities Eq. (|^ we find 


Re. 


f f(u')-5dV 
Jv 

( 

/ S2 
-SireResi^ - 

+ (M • ti) X ti^ • LJ + O(e^Res), 


= -8^Reg( - + 

3 dt 8 


where 



35 


- S°°-S“. (11) 

280 ^ ’ 


So far the hydrodynamic torque acting on the spheroid 
in a general linear flow is given in a form involving the 
scalar product with the angular velocity of the auxil¬ 
iary problem. Since these results are valid for arbitrary 
uj we conclude 


'rh=- 


„ 9 459 n\ „ /3 177 X „ / 13 n\ 

*’'(' - 5 ‘ + 35 o 0 “- + - 350 ' ) ^ ‘0 t. - 8 »(e - jjE ) ((S . t,) X t,) 


177 


13 


-|- 87rReg ( ~ ^ —f g -I- 87reReg ^ ((S°° ■ ti) • 4- (M • ti) x ti^ -|- O(e^Reg) . 


32 


( 12 ) 


This equation constitutes one of the main results of this 
paper. Note that the terms in the first row of the r.h.s 
of Eq. (121 correspond to the creeping-flow limit (RCg = 
0). The terms in the second row account for fluid-inertia 
effects. Since the slip angular velocity scales as 0{e) we 
infer that these effects scale as O(eReg) for any linear 
flow. 


In what follows, an approximate dynamical equation 
for the angular dynamics of the particle is derived. To 
do so, we need to return to Eq. ^ which governs the 
orientational dynamics of the spheroid. The moment-of- 
inertia tensor of the nearly spherical particle reads 


J = 0 ti -I- B^{t2 012 -I- ta (g) ta) 
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where 

= ^(I_4e + 6e2) + 0(e3), 

15 

B' = |(l-3.+ j^“)+0(,=). 

Expanding the time derivative in Eq. (|^ yields 
St^(j • w) = St(J • ^ + X ti)(a; • ti)) . 

Writing this equation in terms of the slip angular velocity, 


and making use of the fact that du)°°/dt = 0 , we are led 
to 


St^fj -u;) = 87rStf4^) 

dtV / V15 dt / 

- 87 reSt(^(a;°° •ti)(a;“ X ti)) 

+ 0(e2St) . 


From Eqs. (§, and ( |13[ ) we obtain 


«=( 


_ „ 459 

5^”^ 350^ 

1 dw, 1 


‘h-{ 


Re, 




3 dt 8 

O^e^St, . 


177 , 

‘^0 ~ l5^‘^ 


) *1 + G - • ti) X ti) + St(-^) 

r 32 


X ti) 


eRes 


'1105 


1 du3, 

15 dt 

ti) • ti )u)°° — (M • ti) X ti 


(14) 


From this equation we compute the particle angular ve¬ 
locity order by order. To this end we insert the ansatz 


u}„ = 


-eu;l'^)-hSt -be^ u}f^ (15) 


,(St). 


,(Re,), 2 „,(.= ) 


-beSta;l®®‘)-beRe,a;i®^®")-bO(e^e2St, e^Re,) , 


into Eq. (14|. To order 0(1) this gives 

= 0. 


This means that to leading order the angular velocity of 
the particle is governed by the vorticity of the unper¬ 
turbed flow (a;l°) = 02 °°). To order 0(e) we find: 

= ti X (S°° • ti) . (16) 

This term is the leading-order term of Jeffery’s angular 
velocity. It combines with the order-e^ term 

a;(®^) = itix(S®®-ti) (17) 

to give Jeffery’s angular velocity Ati x (S°° -ti) to order 
0(e3). To orders 0(St) and O(Res) we find: 

= 0 , = 0 . 


These results suggest that for any linear flow, and up to 
first order, neither particle inertia nor fluid inertia modify 
the orientation of the angular velocity of a sphere as is 
expected from symmetry arguments. Finally, to orders 
O(eSt) and O(eRes) we find 

j^[(0“-ti)x(S®°-ti) 

-bti X (S~ • 0°°-ti) (18) 

-(a;®°-ti) 0 °°-ti] , 


and 


^(.Re.)^__(OOo,t^)x(S-.ti) 

- iti X (S°° • • ti) + ls°° • (ti X (S~ • ti)) 


(19) 


32 

l05 

[(^ 


((S~-ti)-ti) 

3 

35 


3 123 

CXD ^ COO I COO COO 


X t 


1 ■ 


These correction terms account for particle-inertia and 
fluid-inertia effects. The particle-inertia contribution to 
the particle angular velocity was calculated in Ref. [12] 
for a spheroid immersed in a general linear flow. Expand¬ 
ing Eq. ( 8 ) in Ref. [T2| to first order in e results in an 


expression consistent with Eq. (18 1 . 

From Eqs. and we see that both correction 
terms are quadratic in the ambient flow-gradient tensors 
A°°. We also observe that the orientation vector ti oc¬ 
curs twice in each part composing these corrections. This 
implies that the slip angular velocity remains unchanged 
upon replacing ti by — ti (particle inversion symmetry). 

The above results for the angular velocity of the parti¬ 
cle give rise to an effective equation of motion that allows 
to examine the role played by inertial effects on the rota¬ 
tion of a nearly spherical particle in general linear flows. 
To this end we parametrise the orientation of the particle 
using the polar angle 9 and the azimuthal angle (j), see 
Fig. [3 In the laboratory basis the orientation vector 
ti reads 

ti = sin 9 cos (pei + sin 9 sin (j)e 2 + cos 9 63 . (20) 

Eqs. ([3) and (151 - (20) result in a non-linear dynamical 
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system of the form 


in the laboratory basis reads 


^ = f{9,^), e = g{e,4>)- ( 21 ) 


In the following section we discuss the explicit forms of 
this equation for three different linear flows. 


/O 1 0\ 

= 0 0 0 , 
\ 0 0 0 / 


V. EXAMPLES 


Pure shear flow. For a pure shear flow the matrix A°° see Fig. Eq. (21) takes the form 


= -{ecos2(j) — 1) — 


sin 2(j) 

■ f eSt ^ eROs 

\ . 2 n 

37eRes' 

4 

15 35 / 

)sin 0+ 

105 


• sin 29 sin 26 sin 9 cos 9 \ / eSt 
9=e -^^ +^- —- 


LV 15 


37eRes 

105 


eSt lleRCs 


)cos2,f+^+ ^ 


(22a) 

(22b) 


This equation is equivalent to the near-spherical limit of 
Eq. (9) in Ref. [T]. The derivation outlined in the present 
article differs from the calculations in Refs. mm in that 
the method summarised here relies on a basis expansion 
in spherical harmonics, and a joint expansion in the par¬ 
ticle eccentricity e and the shear Reynolds number RCg. 
The calculations in Refs. mm by contrast make use of 
a multipole expansion, valid for arbitrary aspect ratios. 
The fact that the calculations agree lends support to the 
result, to Eq. (221 as well as to Eq. (9) in Ref. [T]. We 
refer the reader to Refs. mm for a further discussion of 
the implications of these results. 

Purely rotational flow. In the case of a purely rota¬ 
tional flow the matrix A°° reads 


/O ~1 0\ 

^oo ^ = 1 0 0 , so that S°° = 0 . 

\ 0 0 0 y 


In this case Jeffery’s slip angular velocity vanishes and 
fluid-inertia effects vanish for small Res. Only particle 
inertia affects the angular velocity of the spheroid, and 
we find: 


</)=! 

9 sin 20 . 
30 


(23a) 

(23b) 


In a rotational flow the evolutions of the Euler-angles are 
thus decoupled. Eq. (23a I indicates that the spheroid 
rotates with the same angular velocity as the fluid. Eq. 
(23b) admits two equilibrium orientations (modulo tt) for 
the angle 9. The first is 0 = 0 (alignment with vorticity), 
and the second is 0 = 7r/2 (the particle rotates in the flow 
plane). Alignment with vorticity is unstable for prolate 
particles (e > 0). For oblate nearly spherical particles we 
find that the stabilities are reversed. 


Purely elongational flow. In the case of a purely elon- 
gational flow, the matrix A°° reads 

(1 0 0 \ 

A°° = = 0 -1 0 , so that = 0 . 

\0 0 oy 

In this case particle-inertia effects vanish so that only 
fluid-inertia effects remain: 

0 = — e sin 2(j) (24a) 

0 = - sin 20 cos 2(5 H- ^ sin 20 . (24b) 

As in the rotational flow the temporal evolution of the 
azimuthal angle (/) decouples from that of the polar angle 
0. Two equilibrium positions are found for the azimuthal 
angle: (p = Q which is stable for a prolate particle (e > 0) 
and unstable for an oblate one (e < 0), and p = 7r/2 
for which the stability is reversed. The equilibrium posi¬ 
tions found for 0 are similar to those obtained for <p, that 
is 0 = 0 and 0 = 7r/2. From the stability analysis of the 
())-dynamics it follows that ecos(2^)/2 remains positive 
for both prolate and oblate particles, so that only the 
polar angle 0 = 7r/2 is found to be stable. As a result, 
prolate and oblate particles orient their axes of symme¬ 
try in the flow plane and finally reach fixed orientation 
repespectively along ±ei (prolate) or along ±e 2 (oblate). 
As Re^ is a small parameter, fluid-inertia effects cannot 
modify these equilibrium positions. They simply speed 
up alignment of prolate particles and slow down align¬ 
ment of oblate particles. 


VI. CONCLUSIONS 

In this paper, an equation of motion has been derived 
for the orientational dynamics of a neutrally buoyant and 
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nearly spherical particle, immersed in a general steady 
linear flow. It would be of interest to extend the results to 
time-dependent linear flows of the form = A°° (t) ■ x. 
In this case the second part of Eq. Q does not hold any 
more and consequently the torque due to unperturbed 
flow does not vanish. Correction terms scaling as O(Res) 


and 0(St) arise. These additional terms are expected to 
render the orientational dynamics of the spheroid more 
complex, but it remains to be seen in detail how these 
terms affect the orientational dynamics of a neutrally 
buoyant particle in an unsteady flow. 
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Appendix: Stokes flow around a spheroid in a 
general linear ambient flow 


In this appendix we describe how to determine the 
Stokes flow around a spheroid immersed in a general lin¬ 
ear ambient flow. To this end, equations Q and (10) 
must be solved. This task is difficult because of the non- 
spherical shape of the particle. This is the reason why 
we have used a perturbation method, briefly described in 
the following. 

The surface of the particle is parametrised by 


r{'d) = 1-1- eh{'d) -\- e^g{'d) 

(see Eq. 0 for the actual definition of the spheroid con¬ 
sidered in this study). We introduce the slip angular 


velocity tensor 

Vt, = Q.-0°° 


where • r = a) x r . 


Performing an expansion to order O(e^) of the fluid ve¬ 
locity in the boundary equation (10) yields 


u' -f (eh{d) 


dW 

dr 


r—1 


(eh(i?) + e^gif})) 




r—1 


dr^ 

= (^1 -b e/i(d) -b 

where = r/|r| is the unit radial vector. We seek a 
solution of Eqs. ([^ and (101 in the form 

u' = Uq -b euj -b e^U 2 -b O(e^) ■ 


Identifying terms of the same order in e results in a set of 
three sub-problems to solve, associated with the bound¬ 
ary conditions 


Ur 


r=l 


u. 


(n, - s~) 

(o, - S“) 
(Qs - S~) 


-h{d) 


dr 


• h{'d)er - 

h{dY 9 ^Uq 

■=i 2 dr"^ r= 


r—1 


r—1 


1 


(25) 

(26) 

(27) 


Lamb solutions. To determine the solutions of Eqs (25) 
to (271, we use general Lamb solutions which exploit the 
fact that any flow solution of the Stokes equations can 
be written in terms of linear combinations of spherical 
harmonics. According to Lamb m (see also Happel & 
Brenner [3]) the general solution for a velocity field that 
vanishes at infinity can be cast in the form 


OO 

u' = ^ [^V X (r X-(n-Hi)) + V$_(„+i) 


n—1 


(n — 2) 2 ^ ^ 

1" V "^-(n-l-l) 


2 n(2n — 1) 
(n-b 1) 


i{2n — 1) 


r Br V- 


T ' -(n+1) 


Here T’ („+!), X-(n-n) and $_(„+i) are spherical har¬ 
monics |13j . These three spherical harmonics involve co¬ 
efficients that remain to be determined using the bound¬ 
ary conditions satished by the fluid velocity on the sur¬ 
face of the particle. Happel and Brenner [3] proposed an 
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efficient method to determine these coefficients that con¬ 
sists in writing the three spherical harmonic functions in 
the form 


is the maximal power of cost'd) or sin(i?) in the bound¬ 
ary condition. The idea here is to fix first the value of 
m of the associated Legendre polynomials by exploiting 
the orthogonality of the trigonometric functions cos{mip) 
and sin(mi^) and then to exploit the orthogonality of the 
Legendre polynomials: 


‘J’-(n-l-l) 


/1\"+1 
(-) [nXn + Y^] 


2(n -I- 1) Vr 


and 


1 

“ n(n + l)VrJ ”■ 

Here X„, and are surface harmonics M- To give 
an example, the surface harmonics read 


Xn =-4„P„(cos(i?)) 


E 


Bnm cos{mip) + Cnm sm{mip) 


(28) 

pr(cos(^)). 


J Pn{x)Pm{x)dx = 


2n + l 

and their associated functions (for a given to): 

f P™(a:)P™da; =_^_ (n -I- to)! 

y_i 2n -I-1 (n — to)! 

As an example, let us show how we determine the co¬ 
efficients An, Bnm and Cnm involved in the surface har¬ 
monic Xn in Eq. (28) from Eq. (29). Consider first the 
case TO = 0. Then: 


Ti = — 
27r 




■ V(i5, </3) )dv? 


Here An, Bnm, and Cnm are constants to be determined, 
and P„ and P™ are the Legendre polynomials and their 
associated functions El- The surface harmonics Yn and 
Zn have expansions similar to Xn but with different con¬ 
stants. The constants are determined by using the rela¬ 
tions 


er-V{d,<p) = J2Xn{^,‘p), (29) 

n—1 


-rX ■Y{d,ip) = Y,Yni^,y^), (30) 

n—1 

and 

OO 

ren-X X V{d,ip) =J2Zni-d,‘p), (31) 

n=l 


where by definition 


r—1 


v{^, . 


and we find for n from 0 to 

An = —-— J sin(il) P„(cos(il)) Tid'd . 
Now consider m > 0. Let 


/■2Tr 


Tim — — 


• V(i?, 1^)^ cos(TO(/3)d:^ , 


and 


72m — 


• V {'d, ip)^ sin(TO(/3)d(p . 


Then we find for n from to to iV: 
2 n -|- 1 (n — to)! 


Bnm — „ / N, I 

2 (n -I- to)! Jo 


and 


Cn.m. — 


2n -I- 1 (n — to)! f 
2 {'n + my.Jo 


sm{'d)Pn{cos{'d)) Timd-d, 


sm{'d)Pn{cos{-d)) T2mdd. 


Now the coefficients in Xn, Yn and Zn are determined as 
follows. Let N be the finite number of harmonics in the 
solution. The value of N does not exceed p -|- 2 where p 


We implemented this method to determine the solutions 
of the three sub-problems associated to the boundary 
conditions (25)-(27) in Maple®. 


















